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We  describe  a  powerful  cross-correlation  technique  for  the  precise  measurement  of  the  proper  motion  of 
tracers  seen  on  successive  images  of  a  time  series  of  solar  granulation.  The  cross  correlation  is  defined  as  a 
function  of  position  in  the  image,  within  a  spatially  localized  apodization  window.fwe  use  a  Gaussian  l"-8"  in 
size).  The  time  average  of  the  spatially  localized  cross  correlation  gives  a  measure  of  the  displacement  that  is 
not  biased  by  atmospheric  seeing.  The  window  size  and  the  seeing  define  the  effective  resolution  of  the  vector 
displacement  determination.  We  use  this  cross-correlation  technique  to  analyze  an  80  minute  run  of  white- 
light  observations  made  at  the  Sacramento  Peak  Vacuum  Tower  Telescope.  Even  though  geometric  distortion 
due  to  atmospheric  seeing  is  instantaneously  at  least  10-20  times  larger  than  the  observed  scale  of  the  large- 
scale  solar  displacements,  100-1000  m/s^^Tits  net  contribution  to  the  80  minute  average  of  proper  motions  is 
RMS  <  20  m/s"1.  The  measured  wclor  displacements  clearly  show  solar  mesogranulation  and  super¬ 
granulation  flows  having  spatial  scales  from  10  to  40^ The  measured  amplitude  of  these  flows  is  significantly 
larger  than  the  RMS  100  nv  s^  noise  which  we  attribute  principally  to  solar  granulation  evolution.  r  •'  : 
Subject  headings :  Sun :  atmospheric  motions/ —  Sun :  granulation  ?  .  r '  • 

,  /  0  it  m  tV  /.->;•  7M  •  • .  l-  . 
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I.  INTRODUCTION 


Proper-motion  measurement  predates  Dopper  spectroscopy  for  measuring  velocities  in  astronomical  objects.  In  astrometry, 
proper-motion  measurement  is  well  defined  when  the  identification  of  a  star  in  two  plates  is  unambiguous  (Stromberg  1946);  star 
catalogs  list  such  measurements  for  most  nearby  stars.  For  the  Sun,  sunspots  have  been  used  as  tracers  to  measure  the  solar  rotation 
rate  (Newton  and  Nunn  1951;  Ward  1 965 ;  Howard,  Gilman,  and  Gilman  1 984). 

Proper-motion  measurements  for  continuous  objects  or  fields  of  stars  is  somewhat  less  well-defined  since  these  objects  evolve  in 
detail  by  changes  in  their  constituent  parts.  Their  proper  motion  can  only  be  given  as  some  form  of  average.  Lindblad  and  Brahde 
(1946)  inferred  the  proper  motions  of  nebular  objects  in  nearby  galaxies.  Simon  (1967),  by  following  the  motions  of  individual 
granules  within  the  time  scale  of  their  lifetimes,  reported  that  granules  tend  to  move  in  a  systematic  way  and  associated  this 
large-scale  movement  with  the  supergranular  convective  flows  (Leighton,  Noyes,  and  Simon  1962;  Simon  and  Leighton  1964). 
Duval  (1980)  applied  cross-correlation  analysis  to  Doppler  velocity  maps  of  the  supergranulation  flows  to  define  the  solar  rotation 
rate. 

Quantified  proper-motion  measurements  of  solar  granulation  have  recently  been  deduced  from  space-based  observations  with 
the  SOUP  experiment  on  Spacelab  2  (Title  et  al.  1987;  November  et  al.  1987).  These  results  were  obtained  using  the  spatially 
localized  cross-correlation  technique  described  by  November  (1986).  The  SOUP  proper-motion  maps  have  revealed  many  new 
features  in  the  quiet  and  active  Sun:  600-1000  ms'1  outflow  from  the  penumbra  and  photosphere  surrounding  a  sunspot,  greatly 
reduced  horizontal  flow  in  a  region  containing  numerous  pores,  quiet  Sun  mesogranulation  and  supergranulation  horizontal 
convective  flow  fields,  and  vortex  flows.  The  proper  motions  are  found  to  be  highly  correlated  with  the  motion  of  magnetic  features 
both  in  direction  and  amplitude  (Simon  et  al.  1988)  and  the  technique  promises  to  give  insight  into  the  foot-point  evolution  and 
twisting  of  solar  magnetic  configurations.  Proper-motion  measurement  of  the  large-scale  horizontal  flow  field  provides  a  unique 
tool  for  the  study  of  solar  convection  that  has  far  greater  potential  than  the  conventional  Doppler  measurement.  The  line-of-sight 
Doppler  velocity  measurement  gives  only  one  component  of  the  velocity  field  whereas  the  proper-motion  measurement  gives  all 
three  when  one  considers  continuity  of  mass  for  a  steady  subsonic  flow  in  a  horizontally  stratified  compressible  fluid  (November  et 
al.  1987). 

Fortunately,  proper-motion  measurement  is  not  restricted  just  to  such  seeing-free  space  observations.  In  this  paper  we  describe  a 
cross-correlation  technique  and  its  application  to  ground-based  observations.  The  atmospheric  geometric  distortion  adds  an  extra 
noise  component  that  is  reduced  by  averaging  methods.  The  basic  technique  is  presented  in  §§  II  and  IV.  The  effect  on  the 
displacement  measurement  caused  by  noise  sources  such  as  geometric  distortion  (atmospheric  seeing  and  film  base  distortion)  and 
blurring  (atmospheric  seeing)  is  determined  theoretically  in  §  III.  A  sample  time  series  of  proper-motion  maps  made  from  ground- 
based  continuum  images  of  solar  granulation  is  analyzed  in  §  V :  The  noise  in  the  measurement  is  treated  using  Fourier  power 
spectrum  analysis  and  determined  as  a  function  of  various  parameters  in  the  measurement  the  time  delay  used  for  the  cross 


1  Operated  by  the  Association  of  Universities  for  Research  in  Astronomy,  Inc.,  under  contract  with  the  National  Science  Foundation.  Partial  support  for  the 
National  Solar  Observatory  is  provided  by  the  United  States  Air  Force  under  a  memorandum  of  understanding  with  the  NSF. 
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correlation,  the  cross-correlation  averaging  time,  and  the  window  size.  This  analysis  permits  us  to  identify  the  noise  sources  and 
their  amplitudes  in  the  data.  Some  of  the  solar  features  that  we  identify  in  our  sample  proper-motion  map  are  discussed  in  §  VI 
along  with  our  conclusions.  However,  this  paper  is  concerned  mainly  with  the  intrinsic  usefulness  and  viability  of  the  technique  in 
observations  that  are  affected  by  atmospheric  seeing.  Thus  most  of  the  observational  results  and  analysis  relevant  to  the  solar 
physics  will  be  discussed  in  subsequent  publications. 

II.  SPATIALLY  LOCALIZED  CROSS  CORRELATION 

We  define  proper  motion  as  the  displacement  that  maximizes  the  spatially  localized  cross  correlation  between  two  images  of  a 
scene  separated  by  a  sampling  time  delay  r  that  is  smaller  than  the  lifetime  of  tracers  in  the  scene.  The  two-dimensional  cross¬ 
correlation  function  is  determined  at  each  “locale”  in  the  scene  by  multiplying  an  intensity  product  with  an  apodizing  window 
W{x).  The  spatially  localized  cross  correlation  C(rf,  x)  is  a  function  of  four  dimensions:  the  two-dimensional  vector  displacement  S 
between  the  images,  and  the  two-dimensional  central  location  x  of  the  window  function.  C{5,  x)  is  defined  in  terms  of  the  intensity 
images  7,(x)  and  7,  +  t(x)  which  sample  the  scene  at  the  2  times  t  and  t  +  r : 

c(s,  x)  =  |  j,(z  -  +  fjmx  -  nd( .  (i) 

The  integral  is  over  the  full  area  of  the  images,  but  effectively  it  is  limited  in  extent  by  the  size  of  the  apodizing  window  function 
W(x).  The  window  size  thus  defines  the  spatial  resolution  of  the  vector  displacement  determination  In  a!!  the  numerical  examples 
used  for  this  work  the  window  was  taken  to  be  a  Gaussian. 

The  spatially  localized  cross  correlation  defined  in  equation  (1)  slides  the  two  intensity  images  symmetrically  in  opposite 
directions.  Equation  (1)  has  the  symmetry  that  the  sign  of  displacement  S  flips  upon  interchange  of  the  images  7,  and  7I  +  t.  Shifting 
the  images  equally  and  oppositely  under  the  window  ensures  that  the  two  image  fields  are  similarly  masked  for  each  spatial 
displacement  5:  The  maximum  of  this  cross  correlation  is  a  most  reasonable  definition  for  the  motion  of  a  tracer  which  undergoes 
evolutionary  changes  between  successive  image  snapshots.  This  formulation  is  not  affected  adversely  by  contrast  variations  that 
occur  over  the  area  of  the  window.  The  images  are  precisely  compared  with  the  same  window  attenuation  in  each  displacement  S.  In 
Appendix  A  we  show  that  our  definition,  equation  (1),  is  simply  related  to  the  conventional  definition  for  the  cross  correlation  in 
which  each  of  the  images  is  independently  masked  to  select  the  area  around  one  spatial  locale. 

Before  forming  the  cross  correlation,  the  images  must  be  spatially  filtered  to  remove  any  large-scale  components.  An  intensity 
gradient  over  the  window  biases  the  cross  correlation  to  give  a  displacement  that  may  not  be  precisely  the  displacement  of  the 
desired  tracer.  Spatial  filtering  to  remove  scales  larger  than  the  window  size  is  a  simple  technique  that  removes  gradients  over  the 
window  function  and  can  make  the  measurement  most  sensitive  to  the  horizontal  transport  of  a  particular  size  of  tracer,  e  g.,  that  of 
the  granulation  field.  We  have  previously  used  Fourier  spatial  filtering  successfully  to  remove  structures  large  compared  to  the 
tracer  (e.g.,  November  et  al.  1987),  but  a  numerically  simpler  and  equally  adequate  filter  is  the  convolution  of  the  intensity  image  /(f) 
with  the  Laplacian  of  a  Gaussian  (LG).  This  yields  the  spatially  filtered  intensity  7(f) :  7(f)  =  i^[/(f)]  =  /(f)  *  [V2  exp  ( -  f  .  f/y2)]. 
This  convolution  can  be  written  as  a  sum  of  separable  components:  7(f)  =  /(f)  *  f(£x)  *  g(fy)  +  /(f)  *  /( fy)  *  g(kx),  where  the 
one-dimensional  /(f)  =  exp[— (f/yj2]  and  g( f)  =  d2/df2[/(f)].  The  direct  numerical  evaluation  of  this  convolution  requires  fewer 
calculations  with  increasing  image  size  than  the  equivalent  Fourier  filter  and  is  easily  implemented  by  direct  integration  on  an  array 
processor.  In  Figure  1  we  show  an  intensity  image  and  its  Laplacian  convolution  for  Gaussian  full  widths  at  half-maximum 
(FWHM)  076,  l'.'O,  174  all  plotted  on  the  same  scale  of  contrast.  The  spatial  filter  having  FWHM  of  170  satisfactorily  removes  all 
scales  except  that  of  granulation  and  is  the  size  that  has  been  used  throughout  this  paper.  The  Laplacian  convolution  also  removes 
power  at  scales  smaller  than  the  diffraction  limit  to  reduce  noise  in  the  proper-motion  measurement.  The  Laplacian  convolution 
does  not  bias  the  displacement  determination;  i.e.,  7(f  +  5)  =  i?[/( f  +  5)]. 

In  order  to  reduce  atmospheric  seeing  effects,  it  is  favorable  to  average  the  cross  correlation  in  time  before  locating  its  maximum 
rather  than  adding  individual  displacement  maps.  The  time-averaged  spatially  localized  cross  correlation  <C(6,  x)>,  is  most  affected 
by  images  in  the  time  series  that  have  the  highest  contrast,  and  we  will  show  in  §  III  that  such  averaging  is  an  optimum  statistical 
procedure  for  adding  low-contrast  contributions  when  the  seeing  is  variable.  An  efficient  numerical  method  to  compute  <C(£,x)>,, 
including  spatial  filtering  by  convolution  with  the  LG,  7(x)  =  i^[/(x)],  is  deduced  by  rewriting  equation  (1)  as  a  spatial  convolution. 
The  spatial  convolution  in  the  multidimensional  x  is  written  as  the  volume  integral :f(x)  *  y(x)  =  |/(f)y(x  -  f)7f.  We  obtain: 

<C<*  x)>,  =  (if  I,(x  -  *  mx)  •  (2) 

Each  temporal  sample  /,(x)  is  spatially  filtered  and  then  multiplied  by  the  spatially  filtered  sample  at  the  fixed  time  delay  t  at  discrete 
spatial  displacements  6.  Typically,  3x3  spatial  displacement  positions  are  used  corresponding  to  dx  =  0,  ±  2  pixels  and  <5y  =  0.  ±  2 
pixels.  (The  smallest  spatial  displacement  increment  that  can  be  easily  represented  is  2  pixels,  corresponding  to  a  shift  of  one  image 
by  one  pixel  in  one  direction  and  a  shift  of  the  second  image  by  one  pixel  in  the  opposite  direction.)  The  image  products  are  averaged 
in  time,  denoted  by  <  ),,  and  the  result  is  smoothed  by  convolution  with  the  sampling  window  IV(x).  We  summarize  schematically 
in  Figure  2  the  procedure  for  calculating  <C(5,  x)),. 

Interpolation  in  the  cross  correlation  <C(5,  x)>,  in  6  is  then  performed  to  define  the  displacement,  A(x),  that  locally  maximizes  the 
cross  correlation  in  each  component  of  S: 
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Fig.  2. — Schematic  representation  of  the  numerical  method  for  the  spatially  localized  cross  correlation  and  displacement  map.  Two  original  granulation  J(x ) 
images  obtained  at  times  t  and  t  +  r  are  spatially  filtered  using  the  convolution  with  an  LG.  The  images  are  shifted  oppositely  and  multiplied  at  nine  spatial  lags  to 
form  nine  cross  product  images  which  are  averaged  in  time  over  the  data  set.  The  average  cross  product  is  smoothed  by  convolution  with  a  Gaussian  window  W  to 
form  the  spatially  localized  cross  correlation  < C(S ,  *)>,.  Using  a  method  of  interpolation  the  maximum  of  the  cross  correlation  is  then  determined  at  each  position  x 
to  give  the  vector  displacement  map  A(jc). 


Usually  it  is  the  case  that  the  image  displacements  we  wish  to  measure  are  much  smaller  than  the  distance  between  pixels  in  our 
digitization  grid.  To  find  an  optimum  interpolation  method  is  not  an  insignificant  task.  We  will  discuss  numerical  procedures  for 
interpolation  in  §  IV  and  in  Appendix  C.  The  spatial  average  of  displacements  over  the  image  field  is  taken  to  be  zero  since  the  exact 
relative  motion  of  the  image  in  the  frame  of  the  object  is  unknown. 

III.  EFFECT  OF  ATMOSPHERIC  SEEING 

The  effect  of  atmospheric  seeing  is  large  in  a  ground-based  observation,  so  it  is  important  that  its  contribution  be  minimized  in 
the  average  of  many  measurements.  The  instantaneous  amplitude  of  the  geometric  distortion  due  to  seeing  may  be  as  large  as  0'.'2  for 
solar  images,  even  in  periods  of  good  seeing  (von  der  Liihe  1985;  November  1986),  and  this  is  at  least  10-20  times  larger  than  the 
amplitude  of  the  flows  we  wish  to  measure.  The  obvious  process  of  temporally  averaging  the  time  series  of  displacement  maps  does 
not  add  the  contributions  from  the  seeing  noise  in  an  optimum  way  when  seeing  conditions  are  variable.  The  local  displacement  is  a 
nonlinear  function  of  the  spatially  localized  cross  correlation  and  is  not  well  defined  when  the  granulation  contrast  becomes  small 
Since  all  the  displacement  measurements  receive  equal  weight  in  the  time-average  displacement  map,  this  procedure  is  not  very 
satisfactory. 

On  the  other  hand,  the  temporal  average  of  the  spatially  localized  cross  correlation  is  affected  in  a  minimum  way  by  the 
poor-seeing  contributions.  When  the  seeing  is  poor  and  the  granulation  contrast  low,  the  cross-correlation  peak  is  broad  and  is  a 
nearly  constant  function  of  displacement.  The  poor-seeing  contribution  to  the  time-averaged  cross  correlation  does  not  bias  its 
centroid  strongly  and  so  gives  a  natural  method  for  weighting  contributions  taken  during  varying  seeing  conditions.  We  consider 
the  effect  that  seeing  has  upon  the  spatially  localized  cross  correlation  function  in  Appendix  B.  We  summarize  our  conclusions  from 
Appendix  B  as  follows:  The  time-averaged  spatially  localized  cross  correlation  including  seeing  is  the  cross  correlation  function 
without  seeing  convolved  with  an  autocorrelation  of  the  time-averaged  seeing  point-spread  function.  The  time-averaged  seeing 
autocorrelation  is  a  symmetric  function  of  displacement  5 ,  and  thus  seeing  tends  not  to  bias  the  cross-correlation  centroid.  A  similar 
cancellation  of  the  seeing  effect  occurs  in  the  displacement  average  when  the  seeing  is  not  variable.  Each  image  is  used  twice  in  the 
average:  the  first  time  an  image  is  used  the  displacement  error  due  to  seeing  adds,  and  the  second  time  it  is  used,  the  displacement 
error  subtracts,  tending  to  offset  the  total  effect.  This  tendency  especially  reduces  the  seeing  effect  in  the  time-averaged  cross 
correlation  because  the  seeing  displacements  tend  to  cancel  even  with  varying  seeing  conditions. 

Nevertheless,  we  thought  it  wise  to  average  the  maximum  number  of  seeing  realizations  in  our  observations.  A  simple  observa¬ 
tional  technique  that  does  this  is  to  take  long  exposures  of  the  granulation  field.  As  long  as  the  image  averaging  time  is  smaller  than 
the  time  delay  t  it  represents  a  viable  sampling.  The  seeing  effect  that  enters  the  spatially  localized  cross  correlation  made  from  a 
small  number  of  long-exposure  solar  images  must  be  the  same  as  that  made  from  a  larger  number  of  short-exposure  images 
spanning  the  same  observation  period.  However,  the  long-exposure  images  have  a  reduced  granulation  contrast,  and  provide  a 
reduced  number  of  independent  contributions  to  the  cross-correlation  average,  which  further  augments  the  loss  of  signal.  In  §  V  we 
demonstrate  the  viability  of  this  procedure  in  a  sample  data  set  taken  with  12  s  exposures  during  a  moderately  good  but  varying 
seeing  condition. 

Another  contributor  to  noise  in  the  displacement  determination  is  geometric  distortion  in  the  film  base.  Film-based  distortion  is 
measured  to  be  less  than  2  pm  (November  1984)  which  corresponds  to  0'.'020  in  the  data  analyzed  in  §  V.  In  a  30  minute  observing 
sequence,  we  take  108  photographic  images  at  the  16.75  s  rate.  This  gives  RMS0.0019  which  is  ~  30  m  s ~ 1  using  a  45  s  time  delay  t 
for  the  cross  correlation. 

IV.  INTERPOLATION  OF  <C(<5,  x)>,  IN  6 

The  cross  correlation  < C(6 ,  x)>,  is  defined  only  coarsely  in  the  variation  in  the  displacement  6  since  the  product  images  from 
equation  (2)  can  be  formed  only  for  displacements  that  are  in  discrete  units  of  the  sampling  interval.  A  precise  definition  for  the 
optimum  displacement  A(x)  that  maximizes  the  cross  correlation  < C(5 ,  x)>„  equation  (3),  requires  interpolation  in  S.  How  closely  we 
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are  able  to  match  interpolated  <C(<5,  x)>,  to  its  actual  functional  form  ultimately  limits  the  precision  of  the  displacement  measure¬ 
ment. 

In  Appendix  C  we  summarize  our  extensive  search  of  interpolation  methods.  We  obtain  this  general  result:  Methods  based  upon 
polynomials  of  degree  less  than  2  always  tend  to  overestimate  the  shift  for  small  displacements,  and  methods  based  upon  poly¬ 
nomials  of  degree  greater  than  2  always  tend  to  underestimate  the  shift.  The  systematic  error  that  occurs  when  the  degree  of 
interpolation  is  other  than  2  is  quite  large,  ~  50%.  Only  quadratic  methods  give  displacements  correct  to  a  precision  of  ~  10%. 

V.  NOISE  ANALYSIS 

We  now  discuss  a  time  series  of  solar  granulation  images  taken  with  12  s  exposure  at  a  16.75  s  rate  in  continuum  (5175.0  A)  at  the 
Sacramento  Peak  70  cm  Vacuum  Tower  Telescope  using  the  Universal  Birefringent  Filter  UBF  (Stauffer  et  al.  1983).  The  observa¬ 
tions  were  made  at  disk  center  under  moderately  good  but  varying  seeing  conditions  for  80  minutes  beginning  at  1645  UT  on  1985 


Fig.  3.— Time-averaged  proper  motion  map  of  solar  granulation.  Vector  displacements  are  shown  superposed  on  a  sample  of  the  solar  granulation  for  an  80 
minute  time  average  of  the  spatially  localized  cross  correlation.  Images  were  spatially  filtered  using  an  LG  with  scale  1,0  arcsec  and  compared  with  the  temporal 
delay  r  =  67  s,  A  spatial  Gaussian  window  of  FWHM  373  was  used  as  indicated  by  the  size  of  the  circle  below  the  figure.  The  scale  of  the  displacements  is  indicated 
by  the  vector  below  the  figure. 


432 


NOVEMBER  AND  SIMON 


Vol.  333 


December  29.  The  images  were  then  digitized  with  the  Sacramento  Peak  flying  laser-spot  Fast  Microdensitometer  (FMS) 
(Arrambide  et  al.  1983)  in  a  spatial  field  64"  x  64"  with  a  Gaussian  spot  of  FWHM  0"33  with  sampling  resolution  of  0"25.  The 
images  were  spatially  filtered  by  convolution  with  an  LG  having  FWHM  F'O.  We  formed  the  spatially  localized  cross  correlation 
using  the  delay  t  =  67  s,  summing  over  the  80  minute  observation,  with  the  Gaussian  window  of  FWHM  373,  following  equation  (2) 
and  the  procedure  shown  in  Figure  2.  We  used  two-dimensional  quadratic  interpolation  in  <C(5,  jc)>,  to  determine  the  displace¬ 
ments  A(x).  Figure  3  is  the  resulting  proper-motion  map.  The  displacements  are  plotted  as  vectors  superposed  on  a  sample 
granulation  image  from  the  sequence.  The  effective  resolution  of  the  displacement  map  defined  by  the  FWHM  of  the  window 
function  is  shown  by  the  circle  in  the  lower  left  of  the  figure.  The  scale  of  amplitude  of  the  displacement  vectors  is  shown  below  the 
figure. 

We  study  the  time  series  of  proper  motion  maps  from  the  data  used  for  Figure  3  in  order  to  understand  the  stochastic  noise  in  the 
measurement.  Proper-motion  maps  were  formed  from  cross  correlations  of  every  third  image  (i.e.,  every  16.75  s  x  3  »:  50  s)  of  the 
observations  using  a  Gaussian  window  of  FWHM  3'.'3  and  two-dimensional  quadratic  interpolation.  A  power  spectrum  was  made 
at  each  point  in  the  spatial  domain  and  in  each  of  its  two  vector  components  by  Fourier  transforming  the  displacement  measure¬ 
ment  in  the  80  minute  time  series  of  96  proper-motion  maps.  The  average  of  these  power  spectra  is  shown  in  Figure  4.  Most  of  the 
power  in  the  spectrum  is  in  the  DC  component  (frequency  0  Hz);  it  is  off  scale  in  the  figure.  For  all  frequencies  other  than  DC,  there 
are  two  identifiable  approximately  spectrally  constant  components:  one  for  frequencies  less  than  4  mHz,  and  one  that  is  evident  at 
higher  frequencies  (above  7  mHz). 

The  standard  deviation  of  a  measurement  is  simply  related  to  measurable  features  of  its  power  spectrum.  We  discuss  this 
relationship  in  Appendix  D  and  demonstrate  that  the  noise  from  an  incoherent  component  of  the  signal  gives  contribution  to  the 
RMS  according  to  its  power.  Assuming  incoherence  between  the  two  spectrally  constant  components  of  Figure  4  we  find  that  the 
low-frequency  component,  with  normalized  power  ~  1.9  x  104(m  s  *)2,  contributes  a  noise  of  ~  140  m  s ' 1  to  the  measurement, 
while  the  high-frequency  component  with  normalized  power  3.6  x  103(m  s1),  contributes  noise  ~60  m  s'1.  The  power  is 
normalized  so  that  it  represents  <s\  according  to  equation  (D3). 

The  low-frequency  component  must  he  due  mainly  to  the  granulation  which  is  known  to  be  evolving  on  a  time  scale  of  ~400  s 
(Title  et  al.  1987).  Spatial  filtering  attenuates  the  5  minute  oscillation  signal  in  its  principal  scale,  and  averaging  over  many  wave 
cycles  attenuates  its  average  signal  in  inverse  proportion  to  the  time,  so  the  5  minute  oscillation  will  contribute  much  less  to  the 
solar  noise  than  the  granulation.  Evolution  of  the  large-scale  flows  on  the  time  scale  of  the  observation  should  also  contribute  to  the 
solar  noise;  however,  no  additional  component  is  evident  in  the  power  spectrum  at  the  lowest  frequencies.  This  indicates  that  the 
large-scale  flows  are  long-lived  compared  to  the  80  minute  time  of  our  observations. 

The  high-frequency  component  appears  to  be  spectrally  white.  We  suspect  that  it  is  largely  caused  by  seeing  effects  and 
photographic  granularity  noise.  However,  a  significant  contribution  may  be  solar  in  origin.  In  Figure  5  we  plot  the  noise  expressed 
in  distance  on  the  Sun  for  the  two  spectral  components  as  a  function  of  the  time  delay  r  used  in  the  cross-correlation  determination. 
For  this  we  reanalyzed  the  data  computing  an  80  minute  time  series  of  48  proper-motion  maps  each  formed  from  averaging  six 
cross-correlation  maps  (100  s).  This  was  done  for  different  delay  times  x  in  units  of  the  16.75  s  sampling  interval.  A  Gaussian  window 


Spatially  Averaged  Temporal  Power  Spectrum 
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Fie.  4. — Spatially  averaged  temporal  power  spectrum  of  displacements 
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Fig.  5.  Displacer. ?nt  noise  as  a  function  of  temporal  delay  The  80  minute  time  scries  of  displacement  maps  was  formed  by  averaging  each  50  s  of  cross 
correlation  maps  using  the  spatial  filter  170  and  the  resolution  window  of  FWHM  5"  for  different  temporal  delay  times  r.  The  RMS  power  in  the  two  identifiable 
spectral  components  from  the  average  temporal  power  spectrum  of  Fig  4  is  shown. 


of  FWHM  5'.'0  was  used  for  all  these  tests.  The  RMS  noise  is  defined  as  the  square  root  of  the  averaged  normalized  power  in  each  of 
the  two  spectral  components,  the  low-frequency  ( <  5  mHz)  and  high-frequency  ( >  S  mHz). 

Detection-noise  contributions  (seeing  and  photographic  granularity)  should  not  change  in  RMS  with  the  temporal  delay  r,  but 
solar-noise  contributions  should  since  the  detection-noise  time  scale  is  short  and  the  solar-noise  time  scale  long  compared  to  the 
delay  time  r.  The  detection-noise  contribution  will  always  be  the  same  between  time  samples  since  the  16.75  s  sampling  time  is  much 
greater  than  the  time  scale  for  the  seeing.  Thus  the  low-frequency  component  is  not  due  principally  to  detection  noise  since  it 
increases  approximately  linearly  with  the  time  delay  t  and  has  zero  intercept  at  r  =  0.  The  high-frequency  component  shows  a 
smaller  increase  with  x  which  then  levels  off  for  r  >  100  s,  and  has  a  positive  intercept  at  x  —  0.  A  conservative  estimate  for  the 
detection-noise  contribution  is  the  noise  at  the  smallest  measured  delay  where  the  solar  features  are  least  changed  by  structure 
evolution;  this  detection-noise  contribution  is  1000  m  on  the  Sun  (less  than  2  miliiarcsec),  or  20  m  s' 1  for  the  time  delay  t  =  50  s. 
The  high-frequency  component  for  time  delays  greater  than  100  s  must  be  mainly  determined  by  those  solar  structures  which  are 
rapidly  evolving,  such  sr  exploding  granules  fTitle  el  al.  1987). 

The  displacement  noise  in  the  low-frequency  component  in  Figure  5  increases  approximately  linearly  with  the  time  delay  t, 
indicating  that  the  solar-noise  contribution  has  an  approximately  constant  signal-to-noise  ratio  independent  of  the  time  delay  t. 
Increasing  the  temporal  baseline  for  the  displacement  measurement  gives  a  proportional  increase  in  the  displacement  signal  and  an 
approximately  proportional  increase  in  the  noise  due  to  the  evolution  of  the  tracer.  In  Figure  6  we  show  the  displacement  noise 
divided  by  the  temporal  delay.  A  time  delay  of  ~r  =  50  s  appears  to  be  the  best  compromise  to  maximize  an  increasing  signal 
relative  to  increasing  noise  due  to  tracer  evolution. 

The  window  size  used  for  the  samples  in  Figures  3  and  4  was  4'.’2  and  in  Figures  5  and  6  was  5'.'0.  For  a  stochastic  process  the  noise 
should  be  inversely  proportional  to  the  window  size  for  spatial  scales  larger  than  either  that  of  the  tracer  or  the  seeing.  We  plot  in 
Figure  7  the  measured  velocity  noise  (solid  line)  in  the  low-frequency  spectral  component  as  a  function  of  the  window  size  for  the  80 
minute  time  series  of  displacement  maps  formed  from  each  100  s  of  data  using  the  time  delay  t  =  67  s.  It  is  clear  that  this  measured 
noise  approaches  asymptotically  the  expected  form  for  a  stochastic  process  ( dashed  line)  as  the  window  size  is  increased. 

We  have  shown  theoretically  in  §  III  that  time  averaging  the  cross  correlation  is  a  more  favorable  procedure  than  time  averaging 
of  displacement  maps.  Figure  8  shows  the  displacement  RMS  power  of  the  low-frequency  component  from  the  temporal  power 
spectra  for  different  cross-correlation  averaging  times  in  the  80  minute  data  set.  For  each  averaging  time  a  time  series  of  displace¬ 
ment  maps  was  formed  with  the  time  delay  t  =  67  s  and  with  the  window  sizeO'.'83  FWHM.  The  noise  decreases  with  the  averaging 
time  and  must  approach  a  constant  value  for  a  large  averaging  time  when  the  nonlinear  effect  of  noise  in  the  displacement 
determination  is  least.  This  figure  indicates  that  little  improvement  can  be  expected  from  displacement  maps  formed  from  cross¬ 
correlation  maps  time  averaged  over  more  than  ~  5  minutes. 

In  earlier  work  we  used  a  technique  for  computing  proper-motion  maps  which  applies  a  sonic  filter  to  white-light  granulation 
images  before  determining  the  displacements  (Title  el  al.  1987;  November  et  al.  1987).  The  sonic  filter  operates  in  the  spatial  and 
temporal  Fourier  transform  domains  of  the  data  to  reduce  to  zero  the  contribution  having  phase  velocity  larger  than  7  km  s  '.By 
this  method  we  eliminated  this  component  which  is  due  principally  to  the  solar  5  minute  oscillation.  We  !■  <  e  applied  this  technique 
to  the  present  data  set  and  find  only  small  differences  between  the  sonic-filtered  map  and  the  time-averaged  proper-motion  map 
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Fir;  6.  Velocity  noise  as  a  function  of  temporal  delay.  The  noise  in  displacement  velocity  is  shown  as  a  function  of  the  temporal  delay  for  the  low-frequency 
component  from  Fig.  5.  The  signal  to  noise  is  best  by  a  small  amount  when  the  temporal  delay  r  =  50  s. 


window  size  f-WHV  :•  >• 

fio  7  Noise  as  a  function  of  window  size  The  RMS  power  from  the  low-frequency  component  in  the  time  senes  of  displacement  maps  is  shown  (sn/uf  /'let  as  a 
function  of  the  resolution  window  size  For  window  size  larger  than  2"  the  noise  approaches  the  form  for  a  stochastic  process  varying  inversely  with  ihe  window  size, 
this  function  is  shown  by  the  dotted  line  For  small  window  sizes  the  sampling  is  both  limited  by  the  scale  of  the  tracer  and  by  the  seeing-disk  size 
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Pig.  8  Noise  as  a  function  of  cross-correlation  averaging  time.  The  low-frequency  contribution  to  the  noise  is  shown  as  a  function  of  the  averaging  time  used  in 
the  spatially  localized  cross  correlation  for  the  displacement  maps  from  the  80  minute  time  series.  The  spatial  filter  was  1.0,  the  temporal  delay  t  =  67  s.  and  Ihe 
resolution  window  size  FWHM  0"83. 


(Fig.  3).  The  temporal  power  spectrum  of  the  displacements  is  similar  to  the  one  we  show  in  Figure  4  in  its  low  frequencies  but 
significantly  attenuated  in  its  high  frequencies.  This  is  expected  as  the  spatial  filtering  of  the  data  by  the  convolution  with  an  LG 
selects  a  specific  range  of  wavenumber.  In  this  range  a  specific  cutoff  frequency  is  defined  by  the  sonic  filter.  The  sonic  filter  does  not 
change  the  actual  noise  but  only  serves  to  extinguish  all  the  temporal  variations  in  the  displacement  map  time  series  above  the  cutoff 
frequency.  In  effect,  sonic  filtering  is  not  different  from  simply  time  averaging  the  cross  correlation,  or  time  averaging  the  displace¬ 
ment  maps.  It  therefore  serves  no  useful  purpose  in  our  analysis  and  will  not  be  used  in  future  investigations. 

VI.  CONCLUSIONS 

We  conclude  that  propter  motions  of  large-scale  solar  features  can  be  determined  reliably  from  long-exposure-time  ground-based 
observations  of  solar  granulation  taken  in  moderately  good  but  varying  seeing  conditions.  The  main  effect  that  atmospheric  seeing 
has  on  the  displacement  measurement  is  to  reduce  the  contrast  in  the  time-averaged  spatially  localized  cross  correlation  without 
biasing  its  central  displacement.  The  reduced  contrast  increases  the  detection  noise  in  the  measurement  but  is  much  smaller  than  the 
noise  introduced  by  granulation  evolution.  In  an  80  minute  average  with  5"  resolution  and  with  the  optimum  delay  of  50  s  there  is 
-  75  m  s  1  noise  due  to  tracer  evolution  compared  to  less  than  20  m  s ' 1  detection  noise. 

We  are  now  obtaining  15  s  exposure  full-disk  observations  of  the  Sun  in  white  light  through  a  6  inch  apjerture  telescopic.  In  times 
of  poor  seeing,  even  though  the  intergranular  lanes  are  not  resolved.  ~  1%  granulation  contrast  is  apparent  everywhere  on  the  disk. 
We  are  optimistic  that  this  experiment  will  be  able  to  provide  full-disk  propier-motion  maps  of  selected  solar  features  or  of  the  full 
disk  on  a  regular  patrol  basis. 

A  large-scale  displacement  map  from  our  data  is  shown  in  Figure  9.  The  proper-motion  map  shown  earlier  (Fig.  3)  lies  within  this 
larger  field,  as  shown  by  the  interior  box.  The  flows  in  the  figure  appear  to  be  long-lived  relative  to  the  80  minute  averaging  time, 
based  upion  the  pxrwer  spectrum  analysis  of  the  time  series.  We  identify  the  coordinates  of  a  number  of  flow  sources  of  amplitude 
500-1000  m  s  1 :  (,x.  >•)  =  (70,  23),  (15,  80),  (55.  120),  (125.  43),  (100,  12),  (100,  108).  We  associate  these  radial  outflows  with  the  40" 
supergranulation  ( Leighton  el  al.  1 962). 

In  Figure  9  there  is  also  much  fine  structure  in  the  horizontal  flows.  This  is  particularly  evident  in  the  divergence  of  the 
displacements  V'  A(x)  which  we  plot  in  Figure  10.  Note  that  positive  and  negative  contours  in  Figure  10  correspond  to  regions  of 
outflow  and  inflow  in  Figure  9,  respectively.  The  contour  interval  is  5  x  10  5  s  l.  and  the  zero  contour  is  not  shown.  The  fine 
structure  has  a  smaller  scale  of  I0"-15"  and  appears  superposed  on  the  larger  flows;  we  associate  this  with  mesogranulation 
(November  el  al.  1981).  Figure  1 1  shows  a  corresponding  contour  plot  of  the  vorticity  of  the  flows.  The  amplitude  of  the  vorticity  is 
roughly  comparable  to  that  of  the  divergence  but  the  features  are  generally  uncorrelated.  The  cyclonic  flows  we  identify  here  must 
be  convective  in  origin.  These  data  are  for  Sun  center;  Coriolis-induced  flows  remain  a  possibility  for  high  latitudes.  The  time  scales 
for  the  flows  is  longer  than  the  80  minute  time  average  for  the  displacement  map.  We  are  now  investigating  the  distribution  of  size 
scales,  the  relation  of  these  flows  to  Doppler  velocity  and  longitudinal  magnetic  field  observations  of  the  Sun.  the  temporal 
evolution  of  the  flow  morphology,  the  fluid  vorticity  and  divergence  relation  for  the  Coriolis  effect  in  the  convection  cells,  and  the 
effective  angular  momentum  transport  over  the  Sun  s  surface  in  these  scales  of  motion.  We  will  discuss  these  issues  in  future 
publications 
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Fig.  9.  —Proper-motion  map  of  a  portion  of  the  solar  image.  It  was  formed  from  the  80  minute  time  average  of  the  spatially  localized  cross  correlation  using  the 
spatial  filter  I  "0,  the  temporal  delay  r  =  67  s.  and  the  resolution  window  3'.'3.  The  scale  of  the  displacements  is  shown  below  the  figure.  This  map  contains  the  smaller 
map  |  Fig.  3)  as  indicated  by  the  enclosed  box.  We  associate  the  largest  amplitude  flows  with  the  supergranular  horizontal  convective  motion. 
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Flti.  10  Divergence  of  flows  The  divergence  is  defined  V  •  AU).  Contours  are  plotted  at  intervals  of  5  x  10  '  s  positive  as  solid  line,  negative  as  dashed  line 
for  the  80 minute  average  from  Fig.  9 (the  zero  contour  is  not  shown).  The  typical  size;.. ale  in  this  map  is  10  15". corresponding  to  the  scale  of  mesogranulalion 
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Fig.  11— Vorticity  of  Flows.  The  vorticily  is  defined  as  £&y(x)/cx  -  ?  &x{x)/?y.  Contours  are  plotted  at  intervals  of  5  x  10  ’  s'  '.positive  as  solid  line,  negative 
as  dashed  line  for  the  80  minute  average  from  Fig.  9  (the  zero  contour  is  not  shown).  The  typical  size  scale  in  the  vorticity  is  also  10”-1 5". 


We  note  in  closing  that  two  other  high-resolution  ground-based  time  series  of  solar  granulation  are  currently  being  analyzed 
using  techniques  similar  to  those  described  herein.  One  uses  the  Swedish  Vacuum  Solar  Telescope  at  La  Palma  (Brandt  el  al.  1988), 
and  the  other  the  Sacramento  Peak  Vacuum  Tower  Telescope  (Title  el  al.  1988).  One  of  their  preliminary  results,  as  in  this  work  and 
in  the  SOUP  observations,  is  that  mesogranulation  is  a  dominant  large-scale  flow  feature  on  the  solar  surface. 

One  of  us  (I.  J.  N.)  initiated  a  study  in  techniques  for  proper-motion  analysis  of  tracers  in  1982  in  collaboration  with  R.  Dunn  and 
R.  Petrov.  The  method  was  first  applied  successfully  in  a  collaboration  between  us  and  T.  Tarbell,  S.  Ferguson,  and  A.  Title  in  the 
analysis  of  Spacelab  2  SOUP  observations  in  1986.  Tarbell  also  demonstrated  to  us  the  usefulness  of  Fourier  techniques  for  noise 
analysis  of  proper-motion  measurements. 


APPENDIX  A 

C(<5,x)  THE  CONVENTIONAL  CROSS  CORRELATION 

Our  definition  for  the  spatially  localized  cross  correlation,  C(S,  x )  from  equation  ( 1 ),  when  multiplied  by  an  attenuation  factor  in  5 , 
W”(dl2),  is  algebraically  equivalent  to  the  conventional  definition  for  the  cross  correlation  of  two  images  in  which  each  is 
independently  masked  to  select  one  locale: 

j  ~  +  +  t)W{x  “  =  j  JW'+M  +  flwa*  -  S)W\X  i-  S)di  .  (Al) 

If  W'(x)  is  a  Gaussian  then  W(x )  and  W"(x)  are  also  Gaussians  having  width  21'2  the  width  of  W\x).  For  the  general  window 
.unction  W'(x)  we  find  W{kl)W"(ki)  =  W'[^(<t,  +  fc2)]IV'[i(*i  —  *2)]  where  W{k)  us  the  Fourier  transform  of  W'(jt) ;  W’{k)  = 
j  W'(x)  exp  (  —  ik'x)dx,  and  kP(/r)  and  W"(k)  are  the  Fourier  transforms  of  the  window  function  W{x)  and  W'(x).  respectively. 
(Herein  we  always  use  the  overbar  to  denote  the  Fourier  transform.)  This  result  is  found  by  substituting  Fourier  integrals  for  the 
window  function  W'{x)  in  equation  (Al)  and  changing  variables  of  integration.  When  the  real  window  function  IV'(j)  is  chosen  to 
have  unit  area,  then  W(<c)  =  W’’(k)  =  W'({k)2.  The  conventional  definition  for  the  cross  correlation  acting  on  images  that  are 
separately  masked,  equation  (Al),  gives  a  slightly  different  measurement  of  the  displacement  than  our  spatially  localized  cross 
correlation  C(d,  x).  The  difference  is  a  displacement  bias  caused  by  the  finite  extent  of  the  window  H/"(<5/2).  Equation  (1)  gives  a  more 
accurate  displacement  measure  because  the  image  fields  are  unbounded  in  x  as  they  are  shifted  in  opposite  directions  under  the  fixed 
sampling  window  W(jt). 
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APPENDIX  B 

THE  EFFECT  OF  SEEING  ON  <C(<5,  x)>, 

We  demonstrate  here  that  seeing  represented  by  a  slowly  spatially  varying  point-spread  function  in  the  object  0,(x)  has  the  effect 
of  convolving  the  time-averaged  spatially  localized  cross  correlation  <C(6,  x)>,  with  a  symmetric  function  in  8.  The  seeing  function  is 
assumed  to  be  uncorrelated  with  the  seeing  function  displaced  in  time  by  the  time  delay  t  for  the  cross  correlation,  and  uncorrelated 
with  the  scene.  We  also  assume  that  the  time-averaged  seeing  function  is  spatially  homogeneous  and  that  the  averaging  time  is  long 
compared  to  the  time  delay  t. 

Let  the  instantaneous  effect  of  seeing  on  the  solar-granulation  intensity  field  0,(x)  be  represented  by  the  following  convolution: 
I^x)  =  J  S,({,  x)0,(x  —  i)d(.  The  instantaneous  seeing  point-spread  function  S({,  x)  is  allowed  to  vary  slowly  in  x  over  the  image. 
The  associativity  property  of  the  normal  convolution  allows  separation  of  the  spatial  filter:  Jif[/t(x)]  =  J  S,({,  x)i^[0,(x  -  {)]d{. 
The  result  is  formally  shown  for  the  slowly  spatial-varying  seeing  point-spread  function  by  writing  i?[/((x)]  as  an  integral,  making  a 
substitution  of  the  variable  of  integration,  and  noting  that  the  spatial  scale  of  the  Gaussian  in  the  integral  (in  this  case  taken  to  be 
the  scale  of  the  solar  granulation)  is  small  compared  to  the  scale  of  spatial  variation  of  the  seeing  point-spread  function.  Thus  we 
have  for  {C(8,  x)>,  from  equation  (2): 

<C(S,  x)>,  =  J  -  <?3)JlV(x  -  tOStf  2,  x)S1  +  t($3,  x)^  d( ,  di2  d( 3  ,  (Bl) 

where  {2,  and  (3  are  variables  of  integration. 

The  time  average  of  a  product  of  incoherent  components  is  the  product  of  the  time  averages  of  each  of  those  components.  Mutual 
incoherence  is  a  reasonable  physical  approximation  between  the  seeing  point-spread  functions  S,  and  S,  +  [  and  the  correlation 
function  if[0,]i^[0,  +  t].  It  is  generally  believed  that  the  seeing  cannot  be  correlated  for  time  scales  longer  than  1  s,  and  in  this 
analysis  the  time  delay  r  is  much  larger,  of  order  1  minute.  Thus  the  seeing  fluctuations  and  the  image  product  must  be  uncorrelated, 
and  we  may  rewrite  equation  (Bl)  as 

<C(8,  x)>,  =  |  (^[°-(^>  -  f  -  +  5  "  <3)  )  W(x  -  {‘)<S(i2 )>'<S(^)>'dfi  dt2di 3 ;  (B2) 

< S(Z )>,  denotes  the  time-averaged  seeing  point-spread  function  which  must  be  the  same  in  the  time  average  <S,>,  and  in  the 
temporally  displaced  time  average  <SI+t>,  in  the  limit  that  the  averaging  time  is  long  compared  to  the  temporal  delay  r.  Also,  we 
assume  that  the  time-averaged  seeing  point-spread  function  does  not  vary  over  the  spatial  domain  of  the  image,  which  must  always 
be  true  in  the  long-time  average. 

The  integral  is  reduced  by  substituting  Fourier  integrals  for  each  of  the  functions:  J^[0,(x)]  =  j  O,  exp  ( ik-x)dk ,  IV(x)  = 
{  W(k)  exp  (ik  •  x)dk,  <£(£)>,  =  f  <$(*)>,  exp  (ik  ■  %)dk.  Of  course,  the  functions  (5, (A),  W(k),  and  <5(*)>,  are  defined  as  the  negative- 
signed-i  Fourier  integrals  of  the  variables  iP[0,(x)],  IT(x),  <S({)>„  respectively;  <S(A)>,  is  the  temporally  averaged  seeing 
modulation-transfer  function.  After  substitution  and  rearrangement  of  the  terms  in  equation  (B2),  ^-function  integrals  occur  w  lich 
simplify  the  integral  and  give  the  following  Fourier  integral  for  (C(8,  x)>,  over  the  two  wavenumber  variables  kt  and  kx : 

<C(8,  x)>,  =  1 0,(\kx  +  k})0,  +  tikx  -  kjy.wj 

<S{±kx  +  *,)>,<$(***  -  *  e“*  ' dkt dkx .  (B. 

It  is  useful  to  write  this  Fourier  transform  as  a  convolution  in  the  spatial  variables  d  and  x  [i.e.,/(£,  x)  *  g(6,  x)  s  J /(e,  {)^(£  -  e, 
x  -  i)dtd f]  : 

<C(6,  x)>,  =  |  (0,(\kx  +  Aa)0,+t(i*x  -  kt)\W(kx)  eikt  *  e^  'dktdk, 

*  |  <S($kx  +  fc4)>,<S(i*x  -  *,)>,  eik‘  *  *  dkt  dkx  .  (B4) 

The  first  integral  in  equation  (B4)  is  the  spatially  localized  cross  correlation  without  seeing;  it  is  convolved  in  (£,  x)  with  the 
seeing-function  integral,  an  integral  that  is  defined  entirely  by  the  time-averaged  seeing  point-spread  function.  The  product 
<s(\kx  +  *,)>,<•?(!*,  -  *  a)>,  in  the  seeing-function  integral  is  invariant  with  the  change  of  sign  of  kt  so  that  the  sine  integral  due  to 
exp  ( ikt  ■  <$)  is  zero.  Thus  the  seeing-function  integral  is  a  symmetric  function  of  8  for  all  kx.  This  proves  that  seeing  has  the  eflect  of 
smoothing  the  spatially  localized  cross  correlation  in  8  by  convolution  with  a  symmetric  function.  The  convolution  cannot  change 
the  centroid  of  the  cross  correlation  in  8  but  may  only  change  its  shape;  i.e.,  broaden  the  cross  correlation. 

In  principle,  seeing  might  introduce  a  location  error  for  the  displacement  since  the  seeing-function  integral  in  equation  (B4)  is  not 
necessarily  symmetric  in  x.  However,  if  <S({)>,  is  a  symmetric  function,  and  real  by  definition,  then  (S(k)),  is  r:al  and  symmetric, 
and  the  seeing-function  integral  is  also  symmetric  in  x.  We  note  that  the  convolution  of  the  spatially  localized  cross  correlation  in  x 
has  the  effect  of  changing  the  shape  of  the  window  function.  Thus  seeing  can  reduce  the  effective  spatial  resolution  of  the 
displacement  map. 

If  the  time-averaged  seeing  point-spread  function  <S({)>,  is  a  Gaussian,  then  the  spatial  and  displacement  convolutions  become 
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separable  in  the  seeing-function  integral  of  equation  (B4).  If  <S({)>,  is  a  Gaussian  then  <S(A)>,  is  also  a  Gaussian  and  the  product 
(S{$kx  +  k6)y,(S($kx  -  Aj)>,  is  the  product  of  a  Gaussian  in  kx  and  a  Gaussian  in  kt.  The  spatially  localized  cross  correlation  from 
equation  (B4)  is  the  spatially  localized  cross  correlation  without  seeing,  convolved  in  S  with  a  Gaussian,  and  then  convolved  in  x 
with  the  product  of  a  Gaussian  and  a  window  function. 

The  result  that  the  time-averaged  seeing  integral  is  a  symmetric  function  of  displacement  S  is  a  direct  consequence  of  the  fact  that 
each  image  is  used  twice  in  forming  the  average,  once  at  time  t  and  again  at  time  t  +  r.  The  first  time  an  image  is  used,  the 
displacement  error  due  to  seeing  adds,  and  the  second  time  the  image  is  used,  the  displacement  error  subtracts,  canceling  its  total 
effect.  The  cancelation  of  the  displacement  error  negates  the  effect  of  geometric  distortion  in  the  time-averaged  cross  correlation  in 
the  limit  that  the  averaging  time  is  long  compared  to  the  time  delay  t  for  the  cross  correlation. 


APPENDIX  C 

INTERPOLATION  IN  C(S,  x) 

The  precision  of  a  displacement  measurement  is  limited  by  the  accuracy  of  the  spatial  interpolation  technique  used.  We  were 
surprised  to  learn,  after  trying  several  methods,  that  this  accuracy  varies  widely  with  the  method,  and  has  a  very  significant  effect  on 
the  result.  Thus  we  undertook  an  extensive  investigation  of  interpolation  techniques,  which  we  describe  here. 

We  base  our  empirical  tests  on  interpolation  techniques  upon  digital  images  made  from  one  photographic  sample  of  solar 
granulation.  The  photograph  was  scanned  with  positional  centers  that  were  displaced  by  different  fractions  of  a  pixel.  The  images 
were  digitized  with  the  Sacramento  Peak  FMS.  The  FMS  has  a  spatial  resolution  of  0.025  in  the  scale  of  the  photographic  sample 
with  positional  accuracy  0''003  mas  (±0.1  pixel).  Two  digital  images  with  centers  displaced  by  0.1  pixel  were  made  from  a  field  of 
granulation  64"  x  64"  using  a  Gaussian  spot  of  FWHM  0'.'5  with  a  sampling  increment  0'.'25.  The  digital  images  were  cleaned  to 
remove  dust  and  other  artifacts  using  a  numerical  method  (November  1988). 

One  of  the  digital  images  was  shifted  by  interpolation  to  a  new  image  center  and  compared  to  the  other  image.  Figure  12  shows 
the  RMS  of  the  difference  of  the  images  as  a  function  of  the  interpolation  shift  for  different  degrees  of  a  polynomial  interpolation. 
We  see  that  the  degree  3  polynomial  (cubic)  gives  by  a  slight  margin  the  minimum  RMS  difference.  For  the  smallest  interpolation 
degrees,  there  is  best  agreement  with  the  offset  digitized  image  when  a  displacement  greater  than  the  true  shift  is  used  for  the 
interpolation.  We  say  that  the  lower  degree  methods  systematically  overestimate  the  actual  shift  by  as  much  as  10%  for  small  shifts. 
The  higher  degree  (greater  than  5)  methods  systematically  underestimate  the  actual  shift  by  as  much  as  40%  for  small  shifts. 

Other  methods  of  interpolation  behave  similarly:  two-dimensional  cubic  splines  through  9x9  points,  one-dimensional  cubic 
convolution  in  4  x  4  points,  and  two-dimensional  quadratic  interpolation  in  3  x  3  points  (six-point  and  nine-point  formulas)  all 
tend  to  overestimate  the  actual  shift  by  10%-25%.  The  interpolation  shift  error  is  found  to  be  fairly  independent  of  the  sample  size 


Fig.  12. — Accuracy  of  images  artificially  shifted  by  five  methods  of  interpolation.  The  normalized  RMS  difference  between  two  images  which  were  digitized  with 
their  centers  offset  by  0.1  pixel  is  plotted  as  a  function  of  the  interpolation  spatial  shift  that  was  applied  to  one  of  the  images  for  various  degrees  of  polynomial 
interpolation.  The  scale  of  the  RMS  was  normalized  by  the  RMS  image  contrast. 
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or  the  sampling  interval  as  long  as  the  solar  granulation  is  properly  resolved.  However,  the  systematic  error  is  noise  sensitive:  if  we 
use  low-pass  spatially  filtered  images,  then  the  higher  order  methods  less  severely  underestimate  the  shift  than  in  the  unfiltered 
examples  shown  in  Figure  12. 

Fourier  interpolation  should  give  an  exact  result  since,  by  the  Fourier  theorem,  a  function  that  is  band-limited  is  exactly 
represented  by  a  finite  number  of  Fourier  sine  and  cosine  terms.  In  practice,  however,  Fourier  interpolation  is  most  sensitive  to 
sampling  effects.  We  find  that  when  the  image  shift  is  near  ^  pixel  its  RMS  variance  with  the  digitized  offset  image  is  several  times 
larger  than  errors  typical  of  the  local  interpolation  methods.  Fourier  interpolation  methods  have  to  be  rejected  as  unsuitable  for 
accurate  displacement  measurements. 

Our  result,  that  interpolation  may  shift  erroneously  an  image  by  a  significant  amount,  warns  of  the  danger  in  using  interpolation 
to  estimate  the  extreme  position  of  a  function.  To  analyze  these  methods  further,  we  digitized  along  the  45  diagonal  the  sample 
granulation  photograph  at  12  evenly  spaced  offset  positions  separated  by  0.1414  pixel.  The  digitized  images  were  numerically 
cleaned  and  spatially  filtered.  The  spatially  localized  cross  correlations  C(S,  x)  was  computed  between  each  offset  digitized  image 
and  the  undisplaced  digitized  image  with  a  Gaussian  window  of  FWHM  3'.'3  in  3  x  3  offsets  (5  of  <5,  =  0,  +  2and<5y  =  0,  +2  pixels. 
A  solution  for  the  extremum  position,  the  displacement  A(x),  is  given  in  the  3  x  3  array  of  C(<5,  x)  at  each  x  by  the  least-squares  fit  in 
a  two-dimensional  quadratic  function  with  6  degrees  of  freedom.  In  Figure  13  we  show  the  spatial  x  average  of  displacements 
normalized  by  the  digitization  offset  as  a  function  of  the  digitization  offset  for  each  offset  digitized  image  compared  to  the 
undisplaced  digitized  image.  The  spatial  average  of  displacements  is  plotted  as  a  point  with  error  bars  at  ±  1  a.  The  accuracy  of  the 
method  is  ~  10%  for  digitization  offsets  less  than  1  pixel  (the  relevant  range  of  offsets  for  this  work).  A  similar  test  made  of  an  image 
field  containing  a  large  sunspot  shows  similar  accuracy  and  noise;  this  demonstrates  the  viability  of  the  method  even  in  a  field 
containing  large  changes  of  structure  and  feature  contrast. 

We  show  in  Figure  14  the  spatially  averaged  displacement  as  a  function  of  the  digitization  offset  for  three  different  methods  of 
interpolation.  The  solid  line  is  for  the  quadratic  interpolation  from  Figure  13.  The  dashed  line  gives  the  estimated  displacement 
using  a  method  of  linear  interpolation  in  3  x  3  points,  while  the  dotted  line  is  the  estimated  displacement  using  a  polynomial  of 
degree  4  (quartic)  fit  through  5  spatial  points  in  each  direction  in  6.  We  have  tried  other  orders  of  polynomials  in  two-dimensional 
and  in  one-dimensional  separable  forms  fit  through  varying  numbers  of  points,  as  well  as  other  methods  including  two-dimensional 
cubic  splines  and  one-dimensional  separable  cubic  convolution. 

In  summary,  our  investigation  gives  this  general  result:  Methods  based  upon  polynomials  of  degree  less  than  2  always  tend  to 
overestimate  the  shift  for  small  displacements,  and  methods  based  upon  polynomials  of  degree  greater  than  2  always  tend  to 
underestimate  the  shift.  The  systematic  error  that  occurs  when  the  degree  of  interpolation  is  other  than  2  is  quite  large,  -  50%.  Only 
quadratic  methods  give  displacements  correct  to  a  precision  of  ~  10%. 


Fig.  13  —  Precision  of  displacement  measurement  using  quadratic  interpolation.  Images  of  a  scene  of  solar  granulation  were  digitized  at  12  positions  in 
increments  of  0.1  pixel  and  displacement  maps  formed  by  comparison  with  the  unshifted  image.  Average  of  displacements  over  the  area  of  each  displacement  map  is 
shown  with  error  bars  representing  ±  I  a.  The  points  are  normalized  by  the  digitization  offset  and  plotted  as  a  function  of  digitization  offset. 
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Fio.  14. — Comparison  of  interpolation  methods.  The  spatial  average  of  displacements  normalized  by  the  digitization  offset  is  shown  as  a  function  of  the 
digitization  offset  between  image  samples  made  of  a  common  scene  of  solar  granulation  using  three  methods  of  interpolation  in  the  spatially  localized  cross 
correlation.  Interpolation  by  two-dimensional  quadratic  ( solid  line )  gives  the  most  precise  measurement.  Interpolation  by  a  two-dimensional  linear  function  ( dashed 
line)  in  the  3  x  3  spatially  localized  cross  correlation  and  interpolation  by  separable  quartic  polynomial  ( dotted  line)  in  the  5  x  5  spatially  localized  cross  correlation 
are  also  shown. 


APPENDIX  D 

NOISE  ESTIMATION  BY  POWER-SPECTRUM  ANALYSIS 

The  standard  formula  for  the  temporal  RMS  standard  deviation  aA  of  the  signal  A(r)  is  al  =  <A(r)2),  —  <A(f)>2,  where  <  >, 
denotes  the  time-averaged  quantity.  This  formula  is  valid  in  the  special  case  when  the  samples  are  independent  and  randomly 
distributed  about  the  average  signal. 

A  more  general  definition  for  the  standard  deviation  that  considers  the  possibility  that  samples  may  not  be  independent  is  given 
from  power  spectral  analysis  of  the  signal  A(t).  Let  A(co)  be  the  Fourier  transform  of  A(t).  Then  by  the  Fourier  power  theorem  we 
have  for  A(a»): 

j*  A(a>)  A(co)*  dco  =  J  A(f)2  df  .  (Dl) 

The  integrals  are  over  the  two  domains  frequency  co  and  time  t.  Let  us  write  A(co)  as  a  superposition  of  mutually  incoherent  signal 
components:  A(co)  =  £  A/co);  the  summation  is  over  the  index  j  which  counts  individual  signal  components.  We  identify  a  special 
DC  component  which  is  the  signal  integral,  A0  s  T <A(f))t,  where  T  is  the  length  of  the  time  domain.  It  is  zero  at  all  frequencies 
u>  i=-  0.  The  DC  component  is  incoherent  with  respect  to  the  fluctuating  components  if  and  only  if  the  fluctuating  components  are 
randomly  distributed.  This  condition  is  usual  for  a  noise  signal  which  fluctuates  about  an  average.  This  assumption  of  incoherence 
between  the  DC  component  and  the  fluctuating  element  A/co)  is  written:  <A0  A/co)>u  =  0.  The  assumption  of  incoherence  among 
all  the  other  individual  elements  obeys  the  same  rule:  <A/co)A*(co)>w  =  0  for  all  j  #  k.  Substituting  the  summation  of  incoherent 
components  for  A(co)  into  equation  (Dl)  gives: 


t  j  AjM  A/co)*  dw  =  r«A(t)2>,  -  <A(f)>2) .  (D2) 

We  have  subtracted  the  squared  DC  component  written  as  a  time  average  from  both  sides  of  equation  (Dl)  and  rewritten  the 
integral  in  time  as  an  average  to  give  equation  (D2).  The  summation  counts  n  incoherent  elements  in  addition  to  the  DC  that  we 
have  identified  for  the  system. 

We  recognize  the  difference  of  averages  on  the  right  side  of  equation  (D2)  which  is  the  usual  formula  for  <j\  for  special  cases. 
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Equation  (D2)  suggests  the  following  as  a  more  general  definition  for  al: 

al  =  Y  i  J  dw  •  (D3) 

Consider  now  the  usual  special  case  that  the  spectrum  contains  two  incoherent  components,  the  DC  component  A0  and  a 
spectrally  white  component  A,(<y).  The  spectrally  white  component  is  constant  in  frequency  for  all  frequencies  up  to  the  Nyquist 
sampling  frequency  av  The  integral  in  frequency  must  be  taken  from  -cuN  to  +coN  and  the  right  side  of  equation  (D3)  becomes  the 
constant  power  |  A,  |2  divided  by  the  number  of  samples,  2 T<uN,  using  the  definition  for  the  Nyquist  frequency.  But  from  the  power 
theorem,  the  power  in  the  white  component  is  simply  the  integrated  signal  for  discrete  samples:  |  A,  |2  =  £  [A(r)  -  <A(t)>,]2,  where 
the  summation  is  over  all  the  samples  in  t.  Thus  equation  (D3)  gives  the  usual  formula  for  the  standard  deviation  in  this  special  case. 

Next  we  consider  another  example,  in  which  the  spectrum  again  contains  two  incoherent  components,  the  DC  and  a  spectrally 
white  component,  but  now  this  white  component  has  a  cutoff  at  frequency  coc  that  is  less  than  the  Nyquist  sampling  frequency.  In 
this  case  equation  (D3)  defines  a\  as  the  power  in  the  white  component,  |  A,  l2,  divided  by  the  reduced  number  of  measurements, 
2Tu>c.  The  usual  formula  for  a\  always  gives  a  lesser  estimate  than  equation  (D3)  in  this  special  case,  and  the  estimate  decreases  in 
value  as  the  sampling  frequency  increases  while  all  other  parameters  remain  fixed.  On  the  other  hand,  a\  from  equation  (D3) 
remains  constant  for  sampling  that  is  greater  than  the  cutoff  frequency.  The  standard  deviation  given  by  equation  (D3)  thus 
correctly  represents  the  fluctuation  in  a  signal  which  is  oversampled.  This  special  case  is  relevant  in  our  discussions  concerning  the 
effect  of  granulation  evolution  on  the  error  in  the  determination  of  displacement,  since  it  represents  an  oversampled  fluctuating 
component. 


REFERENCES 


Arrambide,  M.  R„  Dunn,  R.  B.,  Healy,  A.  W.,  Porter,  R.,  Widener,  A.  L., 
November  ,  L.  J.,  and  Spence,  G.  E.  1983,  in  Proc.  Microdensitometry  Conf. 
Astronomy  (Goddard  Space  Flight  Center),  ed.  D.  Klinglesmith,  NASA  Pub. 
2317,  p.293. 

Brandt,  P.  1988,  Nature,  submitted. 

Duvall,  T.  L.,  Jr.  1980,  Solar  Phys.,  66. 213-221. 

Howard,  R.,  Gilman,  P.  A.,  and  Gilman,  P.  1. 1984,  Ap.  J.,  283,  373. 

Leighton,  R.  B.,  Noyes,  R.  W.,  and  Simon,  G.  W.  1962,  Ap.  J.,  135, 474. 
Lindblad,  B.,  and  Brahde,  R.  1946,  Ap.  J.,  104, 211. 

Newton,  H.  W.,and  Nunn,  M.  L.  1951,  M.N.RA.S.,  111,413. 

November.  L.  j.  1984,  in  Proc.  Sacramento  Peak  Obs..  Small-Scale  Dynamical 
Processes  in  Quiet  Stellar  Atmospheres,  ed.  S.  L.  Keil,  (Sunspot,  NM: 
AURA),  p.  74. 

— - .  1986,  Appl.  Optics,  25, 392. 

- .  1 988,  J.  Opt.  Soc.  Am.,  5, 351. 

November,  L.  J.,  Simon,  G.  W.,  Tarbell,  T.  D.,  Title,  A.  M.  and  Ferguson,  S.  H. 
1987,  in  Proc.  2nd  Workshop  on  Theoretical  Problems  in  High-Resolution 
Solar  Physics,  ed.  G.  Athay  (NASA  Conf.  Pub.  2983),  p.  121. 


November,  L.  J.,  Toomre,  J.,  Gebbie,  K.  B„  and  Simon,  G.  W.  1981,  Ap.  J. 
(Letters),  245,  L123. 

Simon,  G.  W.,  and  Leighton,  R.  1964,  Ap.  J.,  140, 1 1 20. 

Simon,  G.  W.,  1967,  Zs.  Ap.,  65, 345. 

Simon,  G.  W.,  et  al.  1988,  Ap.  J..  327, 967. 

Stauffer,  F.  R.,  Smith,  G.  H.,  Streander,  G.  W„  Wilkins,  L.  M.,  November,  L.  J., 
and  Smartt,  R.  N.  1983,  in  Proc.  Los  Alamos  Conf.  Optics,  ed.  R.  S. 
McDowell  and  S.  C.  Stotlar  (Proc.  SPIE,  290,  p.  77. 

Stromberg,  G.  1946,  Ap.  J.,  104, 12. 

Title,  A.  M„  et  al.  1987,  Adv.  Space  Res.,  68.  No.  8, 253. 

Title,  A.  M.,  Tarbell,  T.  D.,  Shine,  R.  A.,  Topka,  K.  P„  and  Ferguson,  S.  H. 
1988,  Ap.  J.,  submitted. 

von  der  Liihe,  O.  1985,  Lecture  Notes  in  Physics,  Vol.  233,  ed.  R.  Muller 
(Berlin:  Springer- Verlag),  p.  62. 

Ward,  F.  1965,  Ap.  J.,  141,534 


Laurence  J.  November  :  National  Solar  Observatory/Sacramento  Peak  Sunspot,  NM  88349 

George  W.  Simon  :  Air  Force  Geophysics  Laboratory  National  Solar  Observatory /Sacramento  Peak  Sunspot,  NM  88349 


